home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / tmt7.cpp < prev    next >
C/C++ Source or Header  |  1995-01-16  |  3KB  |  107 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5.  
  6. #include "newmat.h"
  7.  
  8.  
  9. /**************************** test program ******************************/
  10.  
  11. void Print(const Matrix& X);
  12. void Print(const UpperTriangularMatrix& X);
  13. void Print(const DiagonalMatrix& X);
  14. void Print(const SymmetricMatrix& X);
  15. void Print(const LowerTriangularMatrix& X);
  16.  
  17. void Clean(Matrix&, Real);
  18.  
  19. void trymat7()
  20. {
  21. //   cout << "\nSeventh test of Matrix package\n";
  22.    Tracer et("Seventh test of Matrix package");
  23.    Exception::PrintTrace(TRUE);
  24.  
  25.    int i,j;
  26.  
  27.  
  28.    DiagonalMatrix D(6);
  29.    UpperTriangularMatrix U(6);
  30.    for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*j-50; D(i,i)=i*i+i-10; }
  31.    LowerTriangularMatrix L=(U*3.0).t();
  32.    SymmetricMatrix S(6);
  33.    for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
  34.    Matrix MD=D; Matrix ML=L; Matrix MU=U;
  35.    Matrix MS=S;
  36.    Matrix M(6,6);
  37.    for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
  38.    {
  39.       Tracer et1("Stage 1");
  40.       Print(Matrix((S-M)-(MS-M)));
  41.       Print(Matrix((-M-S)+(MS+M)));
  42.       Print(Matrix((U-M)-(MU-M)));
  43.    }
  44.    {
  45.       Tracer et1("Stage 2");
  46.       Print(Matrix((L-M)+(M-ML)));
  47.       Print(Matrix((D-M)+(M-MD)));
  48.       Print(Matrix((D-S)+(MS-MD)));
  49.       Print(Matrix((D-L)+(ML-MD)));
  50.    }
  51.  
  52.    { M=MU.t(); }
  53.    LowerTriangularMatrix LY=D.i()*U.t();
  54.    {
  55.       Tracer et1("Stage 3");
  56.       MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
  57.       L=U.t();
  58.       LY=D.i()*L; MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
  59.    }
  60.    {
  61.       Tracer et1("Stage 4");
  62.       UpperTriangularMatrix UT(11);
  63.       for (int i=1;i<=11;i++) for (int j=i;j<=11;j++) UT(i,j)=i*i+j*3;
  64.       GenericMatrix GM; Matrix X;
  65.       UpperBandMatrix UB(11,3); UB.Inject(UT); UT = UB;
  66.       UpperBandMatrix UB2 = UB / 8;
  67.       GM = UB2-UT/8; X = GM; Print(X);
  68.       SymmetricBandMatrix SB(11,4); SB << (UB + UB.t());
  69.       X = SB - UT - UT.t(); Print(X);
  70.       BandMatrix B = UB + UB.t()*2;
  71.       DiagonalMatrix D; D << B;
  72.       X.ReDimension(1,1); X(1,1) = Trace(B)-Sum(D); Print(X);
  73.       X = SB + 5; Matrix X1=X; X = SP(UB,X); Matrix X2 =UB;
  74.       X1 = (X1.AsDiagonal() * X2.AsDiagonal()).AsRow()-X.AsColumn().t();
  75.       Print(X1);
  76.       X1=SB.t(); X2 = B.t(); X = SB.i() * B - X1.i() * X2.t();
  77.       Clean(X,0.00000001); Print(X);
  78.       X = SB.i(); X = X * B - X1.i() * X2.t();
  79.       Clean(X,0.00000001); Print(X);
  80.       D = 1; X = SB.i() * SB - D; Clean(X,0.00000001); Print(X);
  81.       ColumnVector CV(11);
  82.       CV << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
  83.       D << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
  84.       X = CV.AsDiagonal(); X = X-D; Print(X);
  85.       SymmetricBandMatrix SB1(11,7); SB1 = 5; 
  86.       SymmetricBandMatrix SB2 = SB1 + D;
  87.       X.ReDimension(11,11); X=0;
  88.       for (i=1;i<=11;i++) for (j=1;j<=11;j++)
  89.       {
  90.          if (abs(i-j)<=7) X(i,j)=5;
  91.          if (i==j) X(i,j)+=CV(i);
  92.       }
  93.       SymmetricMatrix SM; SM.ReDimension(11);
  94.       SM=SB; SB = SB+SB2; X1 = SM+X-SB; Print(X1);
  95.       SB2=0; X2=SB2; X1=SB; Print(X2);
  96.       for (i=1;i<=11;i++) SB2.Column(i)<<SB.Column(i);
  97.       X1=X1-SB2; Print(X1);
  98.       X = SB; SB2.ReDimension(11,4); SB2 = SB*5; SB2 = SB + SB2;
  99.       X1 = X*6 - SB2; Print(X1);
  100.       X1 = SP(SB,SB2/3); X1=X1-SP(X,X*2); Print(X1);
  101.       X1 = SP(SB2/6,X*2); X1=X1-SP(X*2,X); Print(X1);
  102.    }
  103.  
  104.  
  105. //   cout << "\nEnd of seventh test\n";
  106. }
  107.